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ABSTRACT 

We constrain the isotropic luminosity function (LF) and formation rate of long 7-ray bursts (GRBs) 
by fitting models jointly to both the observed differential peak flux and redshift distributions. We find 
evidence supporting an evolving LF, where the luminosity scales as (1+z) 15 with an optimal S of 1.0±0.2. 
For a single power-law LF, the best slope is ~ -1.57 with an upper luminosity of 10 53 3 erg s _1 , while the 
best slopes for a double power-law LF are approximately —1.6 and —2.6 with a break luminosity of 10 52 7 
erg s _1 . Our finding implies a jet model intermediate between the universal structured e(9) oc 6~ 2 model 
and the quasi-universal Gaussian structured model. For the uniform jet model our result is compatible 
with an angle distribution between 2° and 15°. Our best constrained GRB formation rate histories 
increase from z=0 to z=2 by a factor of ~ 30 and then continue increasing slightly. We connect these 
histories to that of the cosmic star formation history, and compare with observational inferences up to 
z~6. GRBs could be tracing the cosmic rates of both the normal and obscured star formation regimes. 
We estimate a current GRB event rate in the Milky Way of ~ 5 10~ 5 yr" 1 , and compare it with the 
birthrate of massive close WR+BH binaries with orbital periods of hours. The agreement is rather good 
suggesting that these systems could be the progenitors of the long GRBs. 



1. INTRODUCTION 

The finding of the luminosity function (LF) and forma- 
tion rate history (FRH) of long 7-ray bursts (GRBs) is 
a key issue to understand the nature of these events, as 
well as to use them as potential tracers of the star for- 
mation rate (SFR) in the universe. The most direct way 
to infer both the LF and FRH is based on the observa- 
tional luminosity-redshift diagram (LZD) (Schaefer, Deng 
& Band 2001; Lloyd- Ronning, Fryer & Ramirez- Ruiz 2002; 
Yonetoku et al. 2003). However, this strategy is condi- 
tioned by (i) the bias introduced on the sample by the 
redshift inference method, and (ii) the sensitivity limit re- 
lated to the redshift estimate, which in current LZDs is 
> 1 photons enr 2 s _1 . In order to find a reasonable LF, 
this limit has to be pushed down at least by an order of 
magnitude. These goals will be reached by future space 
missions; for example, the Burst Alert Telescope, BAT on- 
board Swift 1 will be 5 times more sensitive than BATSE, 
reaching a limiting flux of 10 -8 erg cm -2 s _1 . Another 
approach, commonly used in the past, is based on fitting 
models with parametric LFs and with a given FRH to the 



observed peak flux GRB distribution, for which the statis- 
tics is healthy (see Stern et al. 2002a, and Guetta, Piran & 
Waxman 2003, and more references therein). Nevertheless, 
the complicated mixing between the model LF and FRH 
introduces a degeneracy among these factors in the flux 
distribution (e.g., Krumhol, Thorsett & Harrison 1998). 
Some early attempts to fix extra constraints on the fits 
by using the information provided by the few GRBs with 
known rcdshifts were done by Sethi & Bhargavi (2001), 
Stern et al. (2002a,b) and Guetta et al. (2003). 

Taking into account the difficulties mentioned above, 
here we propose a new strategy to constrain the GRB LF 
and FRH from the available observations. We will use 
the widest sample up to date for the differential peak-flux, 
P, distribution (the logN-logP diagram, NPD hereafter) 
(Stern et al. 2002b: SAH02), as well as the differential red- 
shift distribution (the N-z diagram, NZD hereafter) from 
220 GRBs inferred from an empirical luminosity-indicator 
relationship (Fenimore & Ramirez-Ruiz 2000: FR00) or 
from the sample of 33 GRBs with known redshifts (e.g., 
see the compilation given in van Putten & Riegenbau 2003: 
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PR03). Then, through Monte-Carlo simulations we will 
optimize the best LF and FRH models by fitting both 
distributions (NPD+NZD) simultaneously by an accurate 
minimum global % 2 criterion. 

After the discovery of afterglows for the long GRBs, it 
was realized that their emitted luminosity should be colli- 
mated. The simplest model for the jet implies a uniform 
energy distribution across the jet and a sharp drop out- 
side Oj, where Oj is the opening angle related to the ob- 
served achromatic break time in the light curves (Rhoads 
1997). For those bursts with known redshift, and there- 
fore with known isotropic luminosity Lj so , L[ so 9j is nearly 
constant (Guetta et al. 2003), although its distribution re- 
sults wider than the one found for Ei so 0j (Frail et al. 2001; 
Bloom, Frail & Kulkarni 2003). A further structured jet 
model proposes a universal non-uniform energy distribu- 
tion within the jet, the observed diversity of break times 
in the light curves being caused by a variation of the view- 
ing angle 9 V (Postnov, Prokhorov & Lipunov 2001; Rossi, 
Lazzati & Rees 2002; Zhang & Meszaros 2002). Similar 
to the uniform jet model, at least for an energy distribu- 
tion per solid angle e(9 v ) oc 9~ 2 , it was found that Li so 6^ 
is also roughly constant. An alternative approach sup- 
poses a Gaussian energy distribution with angle (Zhang & 
Meszaros 2002; Kumar & Granot 2003). Thus, the ori- 
gin of the isotropic LF of GRBs could be due in large part 
to the diversity of opening or viewing angles. Therefore, 
any improvement on the LF knowledge means a progress 
on the understanding and constraining of the jet structure 
(e.g., Guetta et al. 2003; Nakar, Granot & Guetta 2003). 

The collimation of the luminosity is also important at 
the time to calculate the enhancement factor between the 
observed and the true event rates of GRBs. The enhance- 
ment factor for the uniform jet model has been estimated 
previously to be <500 (Frail et al. 2001; PR03), while 
for the structured jet model, values smaller by more than 
an order of magnitude were obtained (Zhang et al. 2003; 
Guetta et al. 2003). Guetta et al. (2003) have also found 
a small value for this factor (75 ± 25) in the case of the 
uniform jet model. Unlike Frail et al., Guetta et al. used 
a weighted average of the angular distribution. 

An ultimate goal of astronomical studies on GRBs is 
the connection of the GRB FRH to the SFRH in the uni- 
verse. Several pieces of evidence show that long GRBs are 
associated to the core collapse of massive stars. Since the 
detected 7-ray fluxes may come from very high-redshifts 
and from dust enshrouded media, GRBs may offer an in- 
teresting way to explore massive star formation in galax- 
ies without the uncertainties of dust extinction, as well 
as in the universe up to very high rcdshifts (e.g., Totani 
1999; Wijers et al. 1998; Lamb & Reichart 2000; Blain & 
Natarajan 2000; Ramirez-Ruiz, Threntan & Blain 2002; 
Choudhury & Srianand 2002; Bromm & Loeb 2002). 

Here we will attempt to constrain the GRB FRH by 
using both the NPD and NZD, and devoting special at- 
tention to the behavior at z>2. We will also make use 
of our results to explore the nature of the GRB progen- 
itors. The association of long GRBs to SNIb/c narrows 
the diversity of potential progenitors to massive helium 
(Wolf-Rayet, WR) stars, which, owing to the angular mo- 
mentum requirement, should be in a close binary system 
with a compact companion. An important task is to esti- 



mate the birthrate of these systems in the Milky Way and 
compare it with the birthrate inferred for the GRBs. 

In §2 we describe the observational data that will be 
used in our study. In §3 the model and strategy applied 
to constrain the LF and FRH from observations are de- 
scribed. In §4, results regarding the best fitting LFs and 
FRHs are presented; evolving LFs are strongly favored, 
and the FRH shape is found to be topologically similar to 
the global SFR history (SFRH). In §5 we discuss the impli- 
cations of an evolving LF on current jet models, while in §6 
the connection between GRB FR and SFR, and some im- 
plications for the progenitors are discussed. Finally, our 
conclusions are given in §7. Throughout this paper we 
use the flat ACDM cosmology with fi m =0.29, fi A =0.71, 
h=0.71. 

2. THE DATA 

Stern et al. (2001,2002a) have processed a sample of 
3255 BATSE triggered and non-triggered GRBs longer 
than 1 s. By using an appropriate efficiency matrix 
they were able to extend the peak flux limit down to 0.1 
photons cm~ 2 s _1 . In order to extend the flux distribu- 
tion in the bright side, the data from Ulysses satellite 
have been used by SAH02 by doing a cross-calibration 
of the joint Ulysses /BATSE events. The final result is 
the differential peak-flux distribution corrected by a joint 
Ulysses/BATSE exposure factor that extends from P=0.1 
to about 300 photons cm~ 2 s _1 . This is the NPD that we 
use here (kindly made available in electronic form to us 
by B.E. Stern). In some previous works, the cumulative 
P distribution has been used. As noticed by the referee, 
random errors in a cumulative distribution propagate in 
an unknown way, therefore for the analysis we pretend to 
do, the differential distribution should be used. 

For the differential redshift distribution we use two sam- 
ples. The first one is a set of 220 BATSE GRBs with a sen- 
sitivity threshold of 1.5 photons cm~ 2 s _1 and with red- 
shifts inferred by using the luminosity- variability empirical 
relation derived by FR00 (kindly made available in elec- 
tronic form to us by E. Ramirez- Ruiz). The second sample 
is a set of 33 GRBs with known z taken from PR03 (exclud- 
ing GRB980425 and adding GRB030429). This sample is 
highly biased by selection effects. We attempt to correct 
some of these effects following Bloom (2003). We use the 
probability of redshift detection due to observability of 
lines in the spectral range and in the presence of night-sky 
lines estimated by Bloom (2003) (who is kindly acknowl- 
edged by making available for us an electronic table with 
the probabilities). The effects related to the distance are 
accounted with another probability which we set equal to 1 
for z< 1 and decreasing as l/(l+z) for larger redshifts (see 
Gou et al. 2003 for the flux detection of GRB afterglows 
located at different redshifts). 

3. THE MODEL 

The differential NPD is modeled by seeding a large num- 
ber of GRBs with a given FR, pgrb (per unit of comov- 
ing volume), and LF, /(L; so )dLi so , and by propagating the 
flux of each source to the present epoch. The rate of GRBs 
observed with peak fluxes between Pi and P2 is (oc NPD): 

N(P 1 ,P 2 )= f dz ^^M f L2 /(Liso )dL iso , (1) 
Jo az 1 + z Ai 
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with Li=Lj so (Pi,z) and L2=Li so (P2,z), while the rate of 
GRBs observed between zi and z 2 (oc NZD) with a limit- 
ing sensitivity P m in is: 



N(zi,z 2 ) = 



22 , dV(z) p GR b(z) 

dz — - 

dz 1 + z 



/(L iso )dL iso , (2) 



where L m =Li so (P m i n ,z), dV(z)/dz is the comoving volume 
at z, the term l/(l+z) takes care of time dilation, and 
L iso (P,z) is obtained by solving the equation: 



(i+-)/(i+z)E E r; s ( E ) dE 



p = 



47Td£(z) 



(3) 



where E„ 



=50 keV and E r 



=300 keV define the sen- 



sitivity band, S(E) is the source differential rest-frame 
photon luminosity, and dL(z) is the luminosity distance 
(see Porciani & Madau 2000). Li so in the rest frame is 

defined as Li so = J 3 1 ° 1 | ) c C ^ kcV ES(E)dE. For the spectrum, we 
use the Band et al. (1993) proposed fit to a sample of 
BATSE spectra, with the parameters taken from averages 
re-derived by Preece et al. (2000) (a = -1, = -2.25, 
Eb = 256 keV). Notice that these values are in the observer 
rest frame. It is assumed that the shape of the Band et al. 
spectrum does not change with z. In the assumption that 
the average redshift of the GRB sample used here is z=l, 
E b would be 511 keV (Porciani & Madau 2000). 

Recently it was discovered a relation between the ob- 
served peak energy in the spectrum (related to Eb) and 
Li so (Yonetoku et al. 2003; see Amati et al. 2002 and 
Lamb, Donaghy & Graziani 2003 for the analogous corre- 
lation between Eb and Ei so ). We will explore both cases, 
a rest Eb constant and a rest 



E b ~ 15 (L iso /10 5U erg s" 1 ) ' 5 keV. 



(4) 



Once defined the cosmological model, the LF, and the 
FRH, then we calculate the NPD and NZD trough Monte- 
Carlo simulations. We use a single power law (SPL) LF, 



/(L iso ) oc L 



7 

ISO ) 



Li < Lj so < L u , 



(5) 



supported by several theoretical models, or a double power 
law (DPL) LF: 



/(Liso) 



T . _ 7i 

-'-'ISO i 



T . ~72 

-'-'ISO ) 



Li < Liso < Lb 



Lb < Liso < L u 



(6) 



suggested by several studies (e.g., Schmidt 1999,2001; 
FR00; Schaefer et al. 2001; SAH02). From a statisti- 
cal analysis of a sample of GRBs with z and Li so inferred 
from empirical luminosity-estimator relationships, Lloyd- 
Ronning et al. (2002) and Yonetoku et al. (2003), have 
suggested that the luminosity should scale with redshift as 
(1+z) 5 , with S wl.4 and 1.9, respectively. 

The GRB FRH per unit of comoving volume is given 



by: 



Pgrb = K/j SF , 



(7) 



where psf is the GRB FRH normalized to a level close to 
the one of cosmic SFRH, and defined as: 



3 e bz 

Psf = »?(z) ■— M Mpc- 3 yr- 1 . 



(8) 



From observations of cosmic SFRH, Porciani & Madau 
(2001) fit a = 22 and b = 3.4. Varying a, this equation 
allows to study the effects of the SFR decline from z=2 to 
0. The function r?(z)=l for z<2, and r](z)=c c(z - 2 '> for z>2, 
allows to study the effects of a growth (decline) of psf for 
z>2 depending on the value of c > (c < 0). Notice that 
the factor K gives information on the GRB progenitor FR 
and has unities of M Q _1 . 

The strategy is to constrain the LF parameters (7,L U ) 
for SPL (keeping Li=10 49 erg s" 1 ), and (71, Lb, 72) for 
DPL (keeping Li=10 49 erg s _1 and L u =10 55 erg s^ 1 ), as 
well as the a, b, c, S, and K parameters. The constraints 
will be obtained by a joint fit of the model predictions to 
the observed NPD and NZD (including their errors), using 
a global x 2 criterion. The total \ 2 1S the sum of %np an d 
Xnz- I n order to verify the quality of fit, tests have been 
made assuming the total x 2 to be a linear combination of 
Xnp an d Xnz: w ith weight coefficients that scale inversely 
to the number of data (bins) in the NPD and NZD samples, 
respectively, in such a way that each sample conditions 
the model with the same strength. The parameter opti- 
mization is based on the non linear Levemberg-Marquardt 
method to find the least x 2 (Press et al. 1988). In spite 
of the large uncertainties, the redshift distribution of 33 
GRBs with known z will be also used as a complementary 
test. 



4. RESULTS 

We have run several models with a SPL or DPL LF 
(identified by S an D, respectively), adopting Eb=511 keV 
or Eb as a function of Li so according to eq. (4) (Yonetoku 
et al. 2003; identified by P and Q, respectively), and either 
without evolution (8 =0) or including 8 among the param- 
eters to optimize (identified by and E, respectively). The 
parameters of each one of these models were optimized by 
a simultaneous fit to the observed NPD and NZD samples, 
as described in §3. Table 1 gives for each model, identified 
as was specified above, the values of the fitted parameters, 
where the same symbols of §3 were used. The uncertainties 
are the conventional standard deviations. The number of 
degrees of freedom (for the simultaneous fit) corresponding 
to each models given in Table 1 is 43 minus the number of 
parameters shown in Table 1 for the corresponding model. 
The information about the quality of the fit for the models 
is given in Table 2. For each model, this Table shows the 
total x 2 , the goodness-of-fit test given by the probability Q 
to find a new total x 2 exceeding the current value, the par- 
tial Xnp and Xnz as weu as the significance levels, P^p an d 
P^f , of the Kolmogorov-Smirnov test between the fitted 
distributions and the NPD and NZD samples, respectively. 
We have introduced the Kolmogorov-Smirnov test just for 
completeness; in fact, for all the models the obtained val- 
ues of the significance levels are not sufficiently high as 
to disprove the null hypothesis that the expected and ob- 
served cumulative distributions are from the same parent 
distribution. Nevertheless, even if these significance levels 
are affected by noise, they reveal a minor quality of the fit 
for SP0 and DP0 models as well as an intermediate quality 
for SQ0 and DQ0 models. The Q probability shows clearly 
this situation. 
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In Fig. 1 we plot the observed and the modeled NPD 
and NZD. The models with SPL and DPL LFs are on top 
and bottom panels, respectively. Dotted lines are for mod- 
els without evolution (0), while solid lines are for models 
with the optimized evolution parameter <5 (E). Thin and 
thick lines identify the cases P and Q, respectively. The 
data are the points with error bars from SAH02 (NPD) 
and from FR00 (NZD). 

GRB luminosity function and evolution. The fit 

of models to the data at the side of the lowest P's in the 
NPD fixes the value of the K factor, while the slope here 
determines the 7 (SPL) or 71 (DPL) LF parameters. The 
quality of the fit can be estimated actually at intermediate 
and high P's. A scale invariant (without cutoff) power law 
LF should provide a (log-log) NPD close to a straight line 
because the power index keeps memory in the NPD, while 
the FRH does not influence on it (e.g., Loredo & Wasser- 
man 1998). Any break to scale invariance introduces some 
curvature besides of the volume evolution effect. The im- 
print of the FRH on the NPD becomes more and more 
significant as the LF diverges from scale invariance. On 
the other hand, if the LF evolves, then as S increases, the 
LF tends to concentrate again its influence on the NPD, 
and the curvature in the NPD grows. This is evident in 
the NPDs of Fig. 1. Non-evolving LFs show a low NPD 
curvature because L u or Lb for GRBs at different z spread 
on a large P range, while for increasing evolutionary effects 
this spread is reduced and the curvature grows. Finally, 
it is quite obvious that a reduced range in P (e.g., 1-50 
photons cm~ 2 s _1 ) weakens strongly the condition given 
by the NPD, and no robust predictions can be made; the 
range used here is rather large, 0.1-300 photons cm~ 2 s _1 . 
The results regarding non-evolving and evolving LFs are: 

1. For both the SPL and DPL cases, the non-evolving 
LFs fit poorly the data as one may appreciate by the val- 
ues of the x 2 , Q and P KS of Table 2. In the NPD, the 
non-evolving models show an excess of GRBs at high P, 
while in the NZD their fits to the data appear by far to 
be the best ones (Fig. 1). In terms of our comparative 
analysis, the non-evolving LF can be ruled out. 

2. The quality of the fit becomes excellent when an 
evolving LF is adopted. The evolutionary models indicate 
an optimal value of 5 = 1.0 ± 0.2. Care has to be taken, 
however, on the fact that the LF is sensitive also to the 
constraints on the NZD. Therefore, a careful analysis has 
to be done, using NZDs obtained from samples based on 
different z determinations, in order to achieve a more con- 
clusive value for 6. 

3. The application of our method to the observational 
NPD and NPZ samples used here points out to a marginal 
preference of the SPL LF on the DPL LF. For a more con- 
clusive result on this question, an analysis on a complete 
LZD has to be carried out. 

Taking into account these considerations, our prelimi- 
nary conclusion is that an evolving LF seems to be neces- 
sary in order to reproduce both the NPD and NZD. 

GRB formation history. The factor K connects the 
GRB FR (frequency) to the assumed SFR, and it is fixed 
mainly by the side of low P's in the NPD. Notice that in 
Eq. (8), psf was normalized to about 0.2 MQMpc _3 yr _1 
at z=2 for most of the models (see Fig. 2); a different 



normalization would imply an equivalent inverse scaling 
of K. It is important to remark that the GRB FRH takes 
into account the contribution of the overall events above 
Li at each z (according to the given LF), while in the NZD 
the lower limit is determined by the sensitivity threshold 
on P m i n . As seen in Fig. 1, the best models in the NZD 
are those with an evolving LF. Figure 2 shows the SFRHs 
corresponding to the models of Fig. 1, using the same 
line code. The vertical segmented area shows the statis- 
tical uncertainty of the Q evolving models. The result is 
encouraging: the GRB FRHs best constrained from obser- 
vations follow, at least qualitatively, the observed cosmic 
SFRH up to zm6 (see §6 and Fig. 4), showing a steep 
increase from z=0 to z^2 (a factor of ~30) and a a gen- 
tle increasing towards earlier epochs. This result shows 
the potential usefulness of GRBs as tracers of SFRs at 
high rcdshifts. In fact our results suggest some difference 
among the SFRHs inferred from GRBs and rest-frame UV 
luminosity (Fig. 4). This difference is expected if the 
GRB-based SFRH is tracing both normal and obscured 
SF regimes (see §6). 

The effect on the results of assuming whether a con- 
stant Eb (case P) or an Eb depending on L; so according 
to Eq. (4) (Yonetoku et al. 2003; case Q) is not relevant. 
For models with a non-evolving LF the fits improve a little 
when using Eb depending on L; so , but their quality remain 
in any case much worse than that of the models with an 
evolving LF. In our calculations we have assumed, rather 
arbitrarily, the low luminosity cut-off Li=10 49 erg s _1 . 
Actually the fit is invariant with respect to a decrease of 
Li below ~ 10 50 erg s _1 , mainly because the GRBs with 
lower luminosities appear above the sensitivity limit of 0.1 
photons cm~ 2 s _1 in a very reduced volume around the 
observer. The only fitting parameter that changes with Li 
is K, which scales as Lj~ °' 6 for the best models. The high- 
est correlation between the parameters appears between 
Loga and b, for which the Pearson's correlation coefficient 
is 0.85 and ALoga w 0.4 Ab. 

Finally, we have also carried out our analysis using the 
NZD from the de-biased sample of 33 GRBs with known 
rcdshifts (see §2). Some difficulty derives from the fact 
that the sensitivity threshold, P m i n , of this sample is un- 
known. Therefore, we have experimented choosing differ- 
ent thresholds between 1 and 5 photons cm" 2 s _1 , and 
proved that the changes on LF and FRH are not signifi- 
cant. The constrained SFRH and its statistical dispersion 
are shown in Fig. 2 by the shaded area, assuming P m ; n =5 
photons cm~ 2 s _1 , and 6=1. A qualitative agreement 
with results obtained using the NZD from FR00 is seen. 
Thus, from different observational samples we reach the 
same conclusion: the shape of the GRB FRH approximates 
that one of the cosmic SFRH. 

5. LUMINOSITY FUNCTION AND JET MODELS 

As mentioned in the introduction, several pieces of ev- 
idence suggest that the LF of GRBs could be strongly 
related to jet collimation effects. From our analysis, us- 
ing the observed NPD and NZD, we have obtained 7 = 
1.55 ± 0.05 for the SPL model. This result provides some 
evidence against the universal structured jet model with 
e(0) oc 8~ 2 (see also Guetta et al. 2003; Nakar et al. 
2003), which predicts a differential LF ocL~ 2 (Rossi et al. 
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2002; Zhang & Meszaros 2002), and also against the quasi- 
universal Gaussian jet structured model, which predicts a 
differential LF ocL -1 even if some dispersion of parame- 
ter values is allowed (Lloyd-Ronning, Dai & Zhang 2003). 
Notice, however, that the best slope we have obtained is 
intermediate between these two cases. In the framework 
of structured jet models, the low luminosity branch of the 
LF is testing how the power per unit solid angle behaves 
for relatively large angles (i.e. corresponding to small ap- 
parent luminosities). For power law structured jets, with 
e(0) oc 9~ k , the predicted LF is a power law with slope 
7=l+2/k (Zhang & Meszaros 2002). Therefore, 7=1.5 
corresponds to k=4, i.e. a steeper value than the "canon- 
ical" k=2. This behavior might correspond to the wings 
of the jet, and not to the entire jet structure: after all, to 
avoid divergence, e(9) needs to be a much more moderate 
function of 9 for small 9. Note also that a power law de- 
scription of e(9) may be an oversimplification of the real 
case, which could be described instead by a double Gaus- 
sian profile (one for the core of the jet and another for the 
wing, as suggested by the behavior of the GRB 030329 
light-curve, Berger et al. 2003a). Also, the unification of 
GRB and XRF requires e(9) to fall off more steeply than 
9~ 2 at large angles, to not over-produce XRFs (Lamb et 
al. 2003; Zhang et al. 2003). 

For a uniform jet, L iso =2L o /0| (9 < 9j) and L iso =0 
(9 > 9j), being 9j the beam aperture and L the intrinsic 
GRB luminosity (Guetta et al. 2003; see also Frail et al. 
2001; Bloom et al. 2003). The differential isotropic lumi- 
nosity distribution scales with 9j as dNoc 9? dN,-, where 
dN,- is the differential jet angle distribution. Let assume 
dN, oc 9j v d9j and dNoc L; so ~ 7 dL; so . With some elemen- 
tary algebra one obtains that v = 5 — 27 (see also Guetta 
et al. 2003). Thus, in the case of the SPL LF the open- 
ing angle distribution is: dN,=0 (6j < (2L /L U ) 1 ^ 2 ) and 
dNj oc 9j 19±01 d9 j {9j > (2Lo/L u ) 1/2 ), while in the case 
of a DPL LF, the opening angle distribution is: dN, oc 
9] 0±1 A6j (9, < (2L /L b ) 1 /2) and dN, = 9j 1 - 95±om d9 J 
(9, > (2L /L b )V2). 

Now, taking into account the bursts of known redshift, 
peak flux and jet angles (from PR03 and Bloom et al. 
2003), at z^l we estimate an average value of 2L ~0.8 
10 50 erg/s (see Fig. 3) and 9 b = {2L /L h ) 1 l 2 ~ 2.3°. 
The upper limit 9 U of the jet angle distribution should be 
determined by the lower limit of LF. This is not our case 
because of the NPD is not sufficiently extended toward 
low peak fluxes. Then we adopt an upper limit obtained 
by the L; so =2Lo/6'| relation, the Eq. (4), and a low limit 
E b ~ 50 keV on the ground of the BATSE sensitivity. This 
assumption gives 9 U ~ 15°. Based on geometrical consid- 
erations, values for 9 U up to 60° cannot be excluded. For 
the DPL LF case the enhancement factor f e = JdNj / J dN 
gives f e =3/(#6# u ) <~ 280, the uncertainty being at least a 
factor two. 

Our results suggest a mild luminosity evolution, 
LbOc(l+z)' 5 with 5 <~ 1. We have checked that this lu- 
minosity evolution does not contradict the suggestion of 
a universal energy reservoir for GBRs (Frail et al. 2001). 
To this end we have calculated the observed peak lumi- 
nosities of the bursts listed in PR03 using a conversion 
factor between the listed count rate and flux of 1 count 



cm~ 2 s _1 = 8.7 10 -8 erg cm~ 2 s . For the bursts in 
common with Bloom et al. 2003), we could then calcu- 
late the true peak luminosity using the jet opening angles 
listed in that paper. Figure 3 shows the jet opening an- 
gles and the corresponding true energies, as a function of 
redshift, for the 24 bursts listed in Bloom et al. (2003). 
This figure shows also the true luminosity for the 19 GRBs 
listed both in PR03 and Bloom et al. (2003). For illus- 
tration, the dashed line shows the case of an evolution of 
the true luminosity proportional to (1+z). It is clear that 
this evolutionary behavior does not contradict the existing 
data. Note also that the found evolution of the observed 
luminosity could be associated not to the evolution of the 
true luminosity, but to the evolution of the aperture angle 
of the jet. 

6. GRBS AS TRACERS OF THE COSMIC STAR FORMATION 
RATE AND IMPLICATIONS FOR THE PROGENITORS 

The understanding of the SF processes and history is 
a key issue to integrate a global vision of the universe. 
GRBs offer the hope of a deep insight of these processes 
if we will be able to establish a connection between GRBs 
and stellar evolution. Our results, though still uncertain, 
have shown that the GRB FRH resembles qualitatively 
the observed cosmic SFRH. The SFRH traced by the rest- 
frame UV luminosity and corrected by dust obscuration, 
as presented in Giavalisco et al. (2004), is shown in Fig. 
4. In this figure we plot also our best models from Fig. 
2. A significant contribution to the cosmic SFRH could 
come from sources emitting strongly in the rest-frame far 
infrared/submillimetre (e.g., Blain et al. 1999; Dunne, 
Ealcs & Edmunds 2003). These objects are likely dust 
enshrouded SF bursts induced by galaxy collisions that 
follow the major mergers of dark halos at earlier epochs 
(for a recent review on galaxy formation, see e.g., Firmani 
& Avila-Reese 2003). GRBs could be direct tracers of 
the SFRH of these galaxies, although this is still an open 
question (see e.g., Ramirez-Ruiz et al. 2002; Berger et al. 
2003b; Lc Floc'h et al. 2003; Barnard et al. 2003). 

The possibility that GRBs are tracing also the SFRH 
in obscured galaxies might explain the apparent differ- 
ences between the inferred GRB FRH and the observed 
UV-cosmic SFRH in Fig. 4. Nevertheless, our results are 
still uncertain due to the nature of the data used to con- 
struct the NZD, and are not suitable for a quantitative 
exploration on these aspects. In the future, when more 
data will be collected, we hope that the comparison of the 
SFRHs inferred from luminous SFR tracers and from the 
GRB FRs will help to clarify interesting questions related 
to the nature of the GRBs as well as to the SF processes 
in highly-obscured galaxies and in the high-redshift uni- 
verse, where reionization makes difficult the observability 
of typical emission lines associated to SF. 

The value of K= pgrb / Psf (see Table 1), gives the num- 
ber of observable (beaming selected) GRB per unit of gas 
mass transformed into stars. For the preferred models, 
K« 0.5 10~ 7 Mq . This value is uncertain because it de- 
pends on the assumed psf normalization and on the lower 
LF limit Lj. For a present-day SFR of 4 M^yr" 1 for the 
Milky Way, we obtain an event rate of ~ 2 10 -7 yr _1 , 
which should be enhanced by the factor f o =280 (see §5); 
the true event rate is then ~ 5 10~ 5 yr _1 . This rate is 
rather uncertain, at least by a factor 2. It is about fifty 
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times lower than the frequency of SNIb/c in the Milky 
Way, showing that only a small fraction of WR stars ex- 
ploding as SNIb/c give rise to the GRB phenomenon (see 
also Podsiadlowski et al. 2004). 

Within the framework of the popular collapsar model 
for long GRBs, the high angular momentum requirement is 
the main difficulty (e.g., Zhang & Meszaros 2003). In fact, 
the strong mass loss by the star and the core-to-envelope 
(dynamo) magnetic coupling slow down the core rotation 
(Woosley, Zhang & Heger 2002; Spruit 2002; Heger & 
Woosley 2002; Izzard, Ramirez- Ruiz & Tout 2003). In- 
stead, if the pre-supernova is in a binary system, then 
the spin-orbit tidal interactions may increase the rotation 
even to the limit of Kerr black hole (BH) formation after 
SNIb/c explosion, ensuring a centrifugally supported disk. 

Binary systems with a WR star co-rotating with the or- 
bital motion and with a period of hours (which implies a 
massive BH as the companion) could be the progenitor of 
long GRBs (Tutukov & Cherepashuk 2003; Tutukov 2003). 
In fact, the condition for Kerr BH formation is ldkRk ~ c, 
where u)k, and Rk are the angular velocity, and the ra- 
dius of the formed Kerr BH, respectively, and c is the light 
speed. If the pre-supernova WR core have an angular ve- 
locity loc and a radius Rc, then the angular momentum 
conservation leads to wcRp ~ u>kR%- Because of the 
core co- rotation with the orbit motion, u>c ~ 27r/P rbit, 
then a Kerr BH is produced if cP or bit <> 27rRc(Rc / R-if )• 
If R c ~ 0.5 Rq and R K ~10 7 cm, then P orb it <, 7 hr, 
while the separation of the binary system, if the total mass 
~30Mq, is <;6 Rq. Close WR+BH binaries show these 
properties. Such systems exist in nature, and some poten- 
tial candidates are known, for example Cyg X-3 (WN3-7+ 
BH) with P or bit=0.2 days. These GRB progenitor systems 
might be observed as luminous, massive X-ray binaries. 

We can now estimate the FR of massive WR+BH bi- 
naries in the Milky Way. Using the initial star-formation 
function for close binaries (Iben & Tutukov 1984): 

where a; n is the initial system semi-major axis, Mi is the 
initial mass of the primary, and q=M2 /Mi . This function 
has been determined from observations in the solar neigh- 
borhood and implies a roughly uniform distribution of ai n 
for 10<aj n /RQ <10 6 , a mass distribution of the primaries 
close to the Salpeter IMF, and a uniform distribution of q. 
If we adopt Mi ~ 25 Mq, dai n /a; n ~ 0.5, and dq~0.3, the 
estimate of the Kerr BH formation rate in the Milky Way 
is ~ 10~ 4 yr _1 . This rate, whose uncertainty is about a 
factor 3, corresponds to a population of a few close binary 
WR progenitors of Kerr BHs at present in the Galaxy. 
The Kerr BH formation rate derived through the previous 
purely astronomical arguments matches rather well with 
the formation rate we have inferred from GRBs. Thus, 
a self-consistent scenario supports the idea that the GRB 
progenitors should be close binary massive WR stars, pos- 
sibly with a massive BH companion. 

As well as the SF enshrouded in dense dust clouds, sev- 
eral other effects can influence the shape difference be- 
tween the GRB FRH and the cosmic SFRH obtained from 
rest-frame UV luminosity. Metallicity influences the star 
mass loss and evolution, and consequently the energetic 



and collimation of GRB. Even the IMF and binary forma- 
tion function may change with redshift. All these ques- 
tions should be well understood before using the GRB FR 
as a tracer of the cosmic SFR. 

7. CONCLUSIONS 

In this work we have exploited the observational peak- 
flux distribution for more than 3300 GRBs, and the red- 
shift distribution for 220 GRBs inferred from an empirical 
luminosity-indicator relationship in order to constrain 
jointly the LF and FRH of long GRBs. Our analysis al- 
lows us to draw the following conclusions: 

-For single or double power law LFs, the case of non- 
evolving LF fits poorly the data, while evolving LFs fit 
rather well both the NPD and NZD. The best fits are ob- 
tained for an evolution where luminosity scales as (l+z) s , 
being 5 = 1.0 ± 0.2. This evolution increases the prob- 
ability to observe GRBs from very high redshifts. The 
introduction of a dependence of Lj so on Eb (Eq. 4) has lit- 
tle effect on the models, the most important one being the 
improvement of the fit for the models with non-evolving 
LF. The goodness of the fit of models with both the SPL 
and DPL evolving LFs are excellent. The quality of the 
fit for the former is slightly better than for the latter. 

-The best models provide a GRB FRH that approxi- 
mately resembles the observed cosmic SFRH, in particu- 
lar if the potential contribution of the SFR in the obscured 
regime is taken into account. The FRH rises steeply from 
z=0 to z~2, and then increases gently toward higher red- 
shifts. The results are qualitatively similar when using 
a sample of 33 GRBs with known z and adequately de- 
biased from selection effects. More data, in particular in 
the NZD, are necessary to explore in more detail the con- 
nection of GRB FR to the cosmic SFR. 

-For the SPL LF, the best slope from the fits is 7 = 
1.57 ± 0.03. In the understanding that the LF is closely 
related to collimation effects, this result implies an inter- 
mediate case between the universal structured jet model 
with e(9) cx 9~ 2 and the quasi-universal Gaussian struc- 
tured jet model. For the uniform jet model, the jet angle 
distributions corresponding to the best model were calcu- 
lated, giving an indicative range between 2° and 15° at 
z=l. 

-Our best models give a true (after collimation effect 
correction) GRB FR of - 5 10~ 5 yr" 1 in the Milky Way. 
Based on astronomical arguments we have argued that 
such a FR agrees with that of close binary systems consist- 
ing of a WR star and a possible massive BH, with periods 
of hours. These systems are able to generate a massive 
Kerr BH after the SNIb/c explosion of the WR (helium) 
star. The observational counterparts of these potential 
GRB progenitors should be luminous X-ray binaries (e.g., 
Cyg X-3), which are estimated to be only a few at present 
in the Milky Way. 

We thank J.S. Bloom, E. Ramirez-Ruiz and B.E. Stern 
for providing us with data used in this paper, and the 
anonymous referee for helpful comments and questions. 
We also thank G. Malaspina for technical assistance. Sup- 
port for this work was provided by CONACyT grant 
33776-E to V.A. and by COFIN grant 2001022957/004 to 
G.G. 
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Fig. 1. — Top panel: Peak flux differential distribution (NP) with the axis in the top-right part, and rcdshift differential distribution (NZ) 
with the axis in the bottom-left part, both for a SPL LF. Error bars show the NP data from SAH02 and the NZ data according to FR00, 
respectively. Dotted lines are for models without evolution (0), while solid lines are for models (E) with the evolution parameter S optimized. 
Thin and thick lines identify the cases with Et,=511 keV (P) or depending on L; so according to Eq. (4) (case Q), respectively. Dashed 
straight line is the -3/2 uniform distribution (Euclidean) behavior in the NPD. Bottom panel: Same as in top panel but for a DPL LF. 
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Fig. 2. — SFRHs inferred from the GRB FRHs under an opportune normalization. Symbols and models are the same as in Fig. 1. Vertical 
segmented areas show the 1 a uncertainty for the evolutionary models. Shaded areas show the evolving models, including their uncertainties, 
for the NZD data corresponding to 33 GRBs with z known (PR03) and de-biased from selection effects. 
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Fig. 3. — From top to bottom: jet opening angle 9, true emitted energy E tI . U c,50 (in units of 10 50 erg) and true emitted peak luminosity 
^true 50 ( m uru ts °f 10 50 erg s~ 1 ), as a function of the redshift z. The values (and limits, for a total of 24 data points) of the jet opening 
angles and energies are taken from Bloom et al. (2003), while the values from the true luminosities (19 points including limits) have been 
calculated using the peak count rate of the GRB listed in PR03. In the bottom panel the dashed lines is proportional to (1+z) to show that 
the mild evolution we find is compatible with existing data and the suggestion of a universal energy reservoir. 
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Fig. 4. — Comparison between the observed cosmic SFRH traced by the rest-frame UV luminosity and corrected by dust obscuration from 
Giavalisco et al. (2004) (dots with error bars) and the SFRH obtained from GRB FRH opportunely normalized (shaded area). The selected 
models include the uncertainties and correspond to 5=1, and Eq. (4). 
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Table 1 

Best-fit parameters for the LF models 



Model' 


7 


LogL u b 




Loga 


b 


c 


6 


K c 


SPO 
SPE 


1.58±0.04 
1.53±0.03 


4.06±0.07 
3.3±0.1 




1.9±0.2 
1.7±0.3 


2.3±0.4 
2.7±0.6 


0.25±0.04 
0.20±0.03 


0.9±0.1 


10.4 
3.43 


SQO 
SQE 


1.59±0.03 
1.57±0.03 


3.96±0.05 
3.3±0.1 




2.0±0.2 
1.8±0.2 


2.5±0.4 
2.8±0.6 


0.26±0.03 
0.21±0.03 


0.8±0.1 


11.6 
4.93 


Model* 


7i 


LogL b b 


72 


Loga 


b 


c 


S 


K c 


DPO 
DPE 


0.9±0.4 
1.53±0.03 


2.6±0.3 
2.8±0.1 


2.1±0.2 
3.4±0.5 


1.5±0.4 
1.6±0.3 


1.6±0.6 
2.8±0.7 


0.29±0.05 
0.17±0.04 


1.2±0.1 


0.73 
2.85 


DQO 
DQE 


0.8±0.5 
1.57±0.04 


2.4±0.4 
2.6±0.1 


2.1±0.2 
2.5±0.2 


1.8±0.4 
1.7±0.1 


1.8±0.7 
2.6±0.3 


0.28±0.05 
0.15±0.02 


1.1±0.1 


0.78 
4.83 



"The symbol code means: S (SPL), D (DPL), P (E b =511 keV), Q (E b from Eq. 4), (no evolution) and E (evolution) 
^Luminosity unit 10 50 erg s _1 , energy range: 30-10000 keV 

c Number of GRBs per unit mass of gas transforming in stars in units of 10~ 8 MI 1 



Table 2 

Statistical quality parameters of the models 



Model* 


x 2 


Q 


Xnp 


Xnz 


pKS 
r NP 


pKS 


SPO 


84.3 


2c-5 


68.0 


16.3 


0.46 


0.12 


SPE 


35.7 


0.48 


26.9 


8.8 


le-6 


le-6 


SQO 


56.6 


0.02 


42.8 


13.8 


3e-3 


0.02 


SQE 


31.7 


0.67 


24.4 


7.4 


le-6 


3e-3 


DPO 


74.6 


2c-4 


64.3 


10.3 


0.26 


0.02 


DPE 


35.5 


0.44 


28.0 


7.5 


le-6 


5e-5 


DQO 


50.0 


0.06 


39.8 


10.2 


0.02 


4c-3 


DQE 


36.6 


0.40 


28.3 


8.3 


le-6 


5e-5 



"The same symbol code of table 1 



